Brain-Machine Interface Utilizing Interventions to Emphasize Aspects of Neural Variance and Decode Speed and Angle

ABSTRACT

A brain machine interface (BMI) for restoring performance of poorly performing decoders is provided. The BMI has a decoder for decoding neural signals for controlling the brain machine interface. The decoder separates in part neural signals associated with a direction of movement and neural signals associated with a speed of movement of the brain machine interface. The decoder assigns relatively greater weight to the neural signals associated with a direction of movement.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority from U.S. Provisional PatentApplication 61/837,014 filed Jun. 19, 2013, which is incorporated hereinby reference. This application is a Continuation-In-Part of U.S. patentapplication Ser. No. 12/932,070 filed Feb. 17, 2011, which isincorporated herein by reference. U.S. patent application Ser. No.12/932,070 filed Feb. 17, 2011 claims priority from U.S. ProvisionalPatent Application 61/338,460 filed Feb. 18, 2010, which is incorporatedherein by reference.

STATEMENT OF GOVERNMENT SPONSORED SUPPORT

This invention was made with Government support under contract nos. NIHRO1-NS054283, 1DP1OD006409, NS054283 and HD075623 awarded by theNational Institutes of Health (NIH), and under contract nos.N66001-06-C-8005 and N66001-10-C-2010 awarded by the Defense AdvancedResearch Projects Agency (DARPA). The Government has certain rights inthis invention.

FIELD OF THE INVENTION

This invention relates to brain machine interfaces and control ofprosthetic devices using neural signals.

BACKGROUND OF THE INVENTION

Brain-machine interfaces (BMIs) translate action potentials fromcortical neurons into control signals for guiding prosthetic devices,such as computer cursors and robotic limbs, offering disabled patientsgreater interaction with the world. BMIs have recently demonstratedconsiderable promise in proof-of-concept laboratory animal experiments,as well as in human clinical trials. However, two critical barriers tosuccessful translation remain. First, current BMIs move considerablyslower and less accurately than the native arm. Second, they do notsustain performance across hours and days, or across behavioral tasks,without human intervention. The present invention addresses this needfor increased performance and robustness and advances the art of neuralprosthetics.

SUMMARY OF THE INVENTION

In one embodiment of the invention a method is provided for artificialcontrol of a prosthetic device (e.g. a robotic limb, a computer cursor,or the like). A brain machine interface is stored on a computer-readablemedium and executable by a computer. The brain machine interfacecontains a mapping from neural signals to corresponding intentionestimating kinematics of a limb trajectory. The intention estimatingkinematics includes positions and velocities. The neural signals aresignals from cortical neurons related to volitional tasks.

The prosthetic device is controlled by the stored brain machineinterface using recorded neural signals as input to the brain machineinterface and the stored mapping of the brain machine interface todetermine the intention estimating kinematics. The determined intentionestimating kinematics then controls the prosthetic device and results inan executed movement of the prosthetic device.

During the control of the prosthetic device, a modified brain machineinterface is developed by modifying (and storing by the computer) thevectors of the velocities defined in the brain machine interface. Eachof the modifications of the velocity vectors includes changing thedirection of the velocity vector towards an end target of the executedmovement of the prosthetic device. Once on target, the velocity is setto zero. In one example, the velocity vectors are modified and stored atdiscrete intervals over the course of the executed movement of theprosthetic device.

The modified brain machine interface includes a new mapping of theneural signals and the intention estimating kinematics that can now beused to control the prosthetic device using recorded neural brainsignals from a user of the prosthetic device. In one example, theintention estimating kinematics of the original and modified brainmachine interface includes a Kalman filter modeling velocities asintentions and positions as feedback.

In another embodiment of the invention a neural prosthetic is providedhaving a prosthetic device and a controller. The controller isexecutable by a computer and interfaced with the prosthetic device. Thecontroller has encoded (and stored on the computer) a mapping of neuralsignals and kinematics of the prosthetic device. Control kinematics isdetermined by the controller to control the prosthetic device based onneural signals received from a user of the prosthetic device. In oneexample of this device, the controller includes a Kalman filter whichencodes the velocity kinematics as intentions and position kinematics asfeedback.

Another objective of this invention is to significantly increase theperformance of poorly performing brain machine interfaces (BMIs) by (i)emphasizing neural variance associated with directional and speedinformation, (ii) incorporating remembered prior information about theneural data when the BMI did not perform poorly, if available. In oneembodiment, this is accomplished by (i) using principal componentanalysis (PCA) to reduce the dimensionality of the data anddownweighting dimensions, which do not have significant directionaltuning, and (ii) utilizing previously known information regarding thelow-dimensional structure of motor-related cortices. In anotherembodiment, we also decode speed and angular variables.

A computer-implemented algorithm has two components. First, thealgorithm utilizes an intervention on a Kalman filter observation noiseprocess to re-weight certain neural observations. When the neuralobservations are low-dimensional neural trajectories of the neural dataas found via PCA, the intervention can be used to ameliorate directionalbiases in the decoder by downweighting dimensions that have lessinformation about the directional control of the prosthetic device.Second, the algorithm may utilize prior information regarding the statespace that the neural data is projected into. For example, if the spacewas better sampled when more neurons were previously observed, thatspace can be remembered; then, in experiments where far fewer neuronsare observed, these neurons can be noisily projected into the rememberedspace. The computer-implemented algorithm may also output, in additionto position and velocity (as in other embodiments described herein)speed and angular variables.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows according to an exemplary embodiment of the inventionperformance comparison between native arm control 110, the controlalgorithm according to examples of this invention 120, and a currentstate-of-the-art algorithm 130. Representative traces of cursor pathduring center-out-and-back reaches. The dashed lines are the targetdemand boxes for the eight peripheral targets and the central target iswithin the dashed boxes. Subsequent targets alternated between thecenter target and the peripheral target in the sequence indicated by thenumbers shown. The traces shown are continuous for the duration of allsixteen center-out and return reaches.

FIG. 2 shows according to an exemplary embodiment of the inventionillustrations of an online neural control paradigm and control algorithmtraining methodology. The input to the neural control algorithm at timet is an array of spike counts, y_(t), collected from implantedelectrodes. The control algorithm translates the incoming neural datainto a velocity output v_(t).

FIG. 3 shows according to an exemplary embodiment of the invention, theneural control algorithm two-stage training method. Initially, cursorkinematics and concurrent neural activity are collected during handcontrol (i). These kinematics (ii) are regressed against neural activityto generate an initial neural control algorithm. Then, a new set ofcursor kinematics and neural activity are collected from closed loopcontrol with the initial algorithm (iii). The cursor kinematics observedwith this neural control algorithm (shown by arrows 330 in (iv)) are theBMI algorithm estimates of the monkeys intended cursor velocity. Werefine this estimate by rotating the velocities towards the goal (shownin arrows 340 in (iv)—note that only one arrow 340 is indicated forclarity purposes, however any other arrow that is not labeled 330 in ivis considered arrow 340). This set of modified cursor kinematics isregressed against neural activity to generate the modified BMI controlalgorithm, shown running in closed loop in (v).

FIGS. 4A-D show according to an exemplary embodiment of the invention acomparison of the offline and online perspectives for neural prostheticdesign: In (A) neural data y_(t) and arm kinematics x_(t) are collectedto fit the parameters of a neural prosthetic. (B) is a hypothetical plotof parameter setting vs. quality/performance of the resulting system bysuch a fitting procedure, given the assumptions of the offlineperspective, essentially that the maxima for arm reconstructions and BMIcontrol occur with the same parameter settings. In (C) x_(t) is nolonger arm kinematics. Data are collected under closed loop neuralcontrol and x_(t) is the kinematics of the neural cursor. This modelfitting procedure assumes that ideal parameter settings for a neuralprosthetic and arm reconstructions may vary, as indicated inhypothetical plot (D).

FIGS. 5A-B show according to an exemplary embodiment of the inventionmethods steps of collecting an “intention-based” kinematic training set:In (A) the user is engaged in online control with a neural cursor.During each moment in the session, the neural decoder drives the cursorwith a velocity (vector 510). We assume that the user intended thecursor to generate a velocity towards target 530 in that moment, so werotate this vector to generate an estimate of intended velocity (vector520). This new set of kinematics is the modified cursor kinematicstraining set. (B) is an example of this transformation applied tosuccessive cursor updates.

FIG. 6 shows according to an exemplary embodiment of the invention agraphical model representing the assumptions of a Kalman filter. x_(t)and y_(t) are the system state and measurement at time t, respectively.

FIG. 7 shows according to an exemplary embodiment of the invention agraphical representation of the position/velocity Kalman filter used forneural control.

FIG. 8 shows according to an exemplary embodiment of the invention aposition/velocity Kalman filter modeling position feedback throughcausal intervention. In other words, the Kalman filter models velocitiesas intentions and positions as feedback.

FIG. 9 shows according to an exemplary embodiment of the invention atechnique to improve BMI performance after it has degraded significantlydue to loss of neural recordings.

FIGS. 10-11 show according to an exemplary embodiment of the inventionintervention PCA results. The performance of the intervention PCAalgorithm 1110 is superior to the state of the art decoder.

DETAILED DESCRIPTION

In an example of the invention, monkeys were trained to acquire targetswith a cursor controlled by either native arm movement or a neuralcontrol algorithm. FIG. 1 shows representative continuous, uninterruptedcursor movements for three different cursor-control modalities: nativearm control 110, the control method according to an embodiment of theinvention 120, and an example of a present state-of-the-art velocityKalman filter (VKF, 130). Monkeys were required to move the computercursor to a visual target (inside the dashed boxes) and hold the cursorat the target square for 500 ms to successfully complete acenter-out-and-back reaching task and receive a liquid reward. Cursormovements with a controller according to an embodiment of the inventionare straighter, faster, and stop better than those generated with theVKF (quantified data not shown). Furthermore, cursor movements producedusing the controller according to an embodiment of the invention (120)are qualitatively similar to those produced by the native arm (110).

This near native level of performance was achieved by a redefinedcontrol algorithm, as well as its training methodology, using a feedbackcontrol perspective. As shown in FIG. 2, when an animal is engaged in abrain machine interface (BMI) task the prosthetic device (e.g., computercursor as shown, or a robotic arm) constitutes a new physical plant withdynamical properties that differ from those of the native arm. Thesubject controls this new plant by (i) modulating neural signals(y_(t),) which are then (ii) measured and (iii) decoded into a velocityby the control algorithm, which is (iv) used to update the cursor on thescreen, affecting neural signals on subsequent time steps. Viewing theBMI from this feedback-control perspective leads to two new methods. Thefirst method is a modification of BMI training method. The second methodis an alteration to the design of the control algorithm. These are partof the control algorithm and produce the BMI results 120 in FIG. 1.

Previously reported control algorithms incorporate a mathematical modelbuilt to reconstruct native arm movement kinematics, implicitly assumingthat neural control of the arm and prosthetic device are equivalent.However, embodiments of this invention take a different approach andincorporate the fact that the subject's control strategy may actuallyvary when asked to operate different plants with different dynamics.FIG. 3 outlines the two-stage BMI optimization procedure according toexemplary embodiments of the invention. In the first stage, the monkeycontrols the cursor using his arm (i in FIG. 3). An initial model is fitusing arm trajectories and simultaneously recorded neural signals (ii inFIG. 3). The monkey then controls the cursor using the BMI with thisinitial model (iii in FIG. 3). This model's cursor trajectories are thenrefined by incorporating knowledge of the monkey's intention (iv in FIG.3), leading to a final model (modified BMI) that is used for theremainder of the experiment (v in FIG. 3).

The first method according to an embodiment of the invention is to fitthe BMI against a novel set of kinematics. Ideally, instead of fittingneural activity against the kinematics of native arm movements, we fitto the monkey's intended cursor kinematics during closed loop neuralcontrol. Since one would lack explicit access to the monkey'sintentions, we make the assumption that the monkey always wishes to movetowards the target 310, as this is how the native arm moves and is alsothe best strategy for rapidly acquiring rewards. We employ thisassumption to estimate the monkey's intention from cursor movementscollected during the initial control session (iii in FIG. 3).Specifically, as shown by iv in FIG. 3, the original BMI cursor movementvelocity vectors (320) are rotated towards the target (310). When thecursor is on target, we assume an intention of zero velocity. Thesemodified cursor kinematics and the corresponding neural data are used tofit the modified BMI, thereby arriving at the final control algorithm.

It is important to note that this transformation is applied only to thetraining data: during online control the BMI has no prior knowledge ofthe subjects task or goal, or the placement of the targets in theworkspace. On subsequent BMI sessions, step-one (arm control) can beskipped (i.e. ii in FIG. 3), and the monkey can begin at step-two(initial BMI control) (iii in FIG. 3) with model parameters from forexample a previous day. Although the initial quality of BMI cursorcontrol (arrows 330 in iv in FIG. 3) is often lower on the second daywhen starting directly with step-two, performance is restored when themodel is optimized with the modified cursor kinematics (340, iv in FIG.3—only one arrow 340 is indicated for clarity purposes, however anyother arrow that is not labeled 330 is considered arrow 340). Thisapproach was employed for up to five consecutive days, and performanceis on par with that shown in FIG. 1 (110, 120).

The second method attempts to capture the observation that neuralactivity is correlated with not only the velocity of the cursor but alsoits position. Most existing BMIs model a relationship between neuralfiring and either velocity or position. Human clinical trials have shownthat BMIs modeling velocity have higher performance. However, if thecontrol algorithm only models the velocity relationship, then a positionbased increase in firing will result in an increase or decrease incursor velocity, which is an undesired effect. To mitigate this effect,in one embodiment of this invention, we explicitly model velocity as theusers intention, and model position as a signal that is fed back as theuser controls the cursor and sees its behavior. This feedback-modifiedKalman filter allows the user to control the velocity of the cursor withthe measured neural signals while taking into account that this signalis influenced (or corrupted) by the current position of the cursor.Since the BMI knows the position of the cursor it is not estimated fromthe measured neural signals. Instead, the expected contribution ofposition to neural activity is removed. This results in a lessposition-dependent neural data stream, which in turn enables moreaccurate estimation of the intended cursor velocity as evidenced byhigher performance (120 in FIG. 1).

Long-duration, continuous, high-performance operation is central to thesuccessful translation of BMIs to human patients. To this end, and infurther embodiments of this invention, we made three specific designchoices for the BMI, in addition to the two methods described supra.First, we did not employ spike sorting, the process of classifying eachspike based on the shape of its action potential. The goal of thisclassification is to separate single channels composed of spikes frommany neurons into multiple channels to access the spiking activity ofindividual neurons. This is a standard practice in neuroscience and alsoin BMI-design, and can yield the highest information per electrode, butit adds the challenge of tracking each sorted action potential due tochanges in amplitude over time. To avoid these signal instabilities,instead of spike sorting, we counted the number of threshold crossingsper electrode within a time window. Second, while most BMI designs admitneural data from one or more neural integration time windows (e.g., 100ms time bin; multiple 50 ms time bins in linear filter with history asfar back as 1 sec) we elected to use a single time bin of 50 ms. Thiswas determined through a series of experiments with humans using anonline prosthetic simulator (OPS), indicating that shorter time bins arepreferable due to their lower latency and therefore reduced closed-loopfeedback time. Control experiments with our monkeys performing BMI tasksare consistent with these OPS results. Finally, the number of highlydistinguishable single neurons on an electrode array tends to decreaseover time, however multiunit activity can often have BMI relevant tuningThe results shown in FIG. 1 were acquired from arrays several months toyears post implantation. By employing these somewhat older arrays, whichhad relatively few clearly distinguishable single units, we were able toconfirm that threshold-crossing based multi-unit activity together withthe control algorithm according to embodiments described herein is ableto provide high performance for months and years post arrayimplantation.

These levels of sustained performance and robustness demonstrate thepotential to provide functional restoration for patients with a limitedability to move and act upon the world. For patients suffering fromneurological injury and disease, although descending pathways arecompromised, motor cortex may be intact, permitting the use of thisclass of technology. In recent years, a number of brain interfacetechnologies using a variety of signal sources, such as theintra-cortical technology described here, electroencephalography (EEG),and electrocorticography (ECoG), have been developed. The BMI communitycontinues to create options for disabled individuals and assess relativerisk and benefit. The performance and robustness advances demonstratedherein, as well as the system design methodology, should help increasethe clinical viability of intra-cortical based BMIs.

Experimental Methods

In an exemplary embodiment, experiments were conducted with adult malerhesus macaques, implanted with 96 electrode Utah arrays (BlackrockMicrosystems Inc., Salt Lake City, Utah) using standard neurosurgicaltechniques. Monkeys were implanted e.g. 19-27 months or 4-8 months priorto the study. Electrode arrays were implanted in a region spanning thearm representation of the dorsal aspect of premotor cortex (PMd) andprimary motor cortex (M1), as estimated visually from local anatomicallandmarks.

The macaques were trained to make point-to-point reaches in a 2D planewith a virtual cursor controlled by movements of the contralateral armor by the output of a neural decoder. The virtual cursor and targetswere presented in an immersive 3D environment. Hand position data weremeasured with an infrared reflective bead tracking system. Spike countsused in all studies were collected by applying a single negativethreshold, set to 4.5×root mean square of the spike band of each neuralchannel—no spike sorting was used. Behavioral control and decoding wererun on separate PCs with custom code running on the PC platform,communication latencies between these two systems was 3 ms. This systemenabled millisecond timing precision in all computations. Neural datawas initially processed by a recording system and was available to thebehavioral control system within 5±1 ms. Visual presentation wasprovided via two LCD monitors with refresh rates at 120 Hz, yieldingframe updates within 11±3 ms. The monitors were used to form aWheatstone stereograph using two mirrors to fuse the two displays into asingle stereo percept.

Kalman Filter

Many control algorithms, or continuous decoding methods, have beenstudied for neural prosthetics applications such as the Kalman filter.The basic intended application of the Kalman filter is to track thestate of a dynamical system throughout time using noisy measurements.Although we have a model of how dynamics evolve through time, theunderlying system may not be deterministic. If we know the state of thesystem perfectly at time t, our dynamical model only gives us anestimate of the system state at time t+1. We can use the measurements(or observations) of the system to refine our estimate, and the Kalmanfilter provides the method by which these sources of information arefused over time. The filter can be presented from a dynamical Bayesiannetwork (DBN) perspective, and is considered to be one of the simplestDBNs.

A graphical model of the basic DBN representation of the Kalman filteris shown in FIG. 6. For neural prosthetic applications, the system statevector x_(t) is commonly used to represent the kinematic state. Inembodiments of this invention, the state vector represents position andvelocity of the cursor:

(x_(t)=[pos_(t) ^(vert),pos_(t) ^(horiz),vel_(t) ^(vert),vel_(t)^(horiz),1]^(T))

The constant 1 is added to the vector to allow observations to have afixed offset (i.e., baseline firing rate). y_(t) is the measured neuralsignal, which is binned spike counts. The choice of bin width can affectthe quality of prosthetic control: assuming local stationarity, long binwidths can provide a more accurate picture of neural state but withpoorer time resolution. Thus, there is an implicit tradeoff between howquickly the prosthetic can change state and how accurately those statesare estimated. Typical bin widths used in studies range from 10 ms to300 ms. Through online study, we find that shorter bin widths result inbetter performance. The results discussed in this study use 25 or 50 msbin widths.

When applying the standard Kalman filter, the system is modeled withlinear dynamics, a linear relationship between kinematic state andneural observations, and Gaussian distributed noise and uncertainty:

x _(t) =Ax _(t−1) +w _(t)  (1)

y _(t) =Cx _(t) +q _(t)  (2)

where A

R^(p×p) and C∈

represent the state and observation matrices, and w and q are additive,Gaussian noise sources, defined as w_(t)˜N(0, W) and q_(t)˜N(0, Q). A isthe linear transformation from previous kinematic state, or dynamics,and C is mapping from kinematic state to expected observation. Thisformulation allows for very fast inference (decoding) of kinematics fromneural activity and the parameters θ={A, C, W, Q} can be quickly learnedfrom training data. The observation model of the Kalman filter, C and Q,is fit by regressing neural activity on observed arm kinematics:

$\begin{matrix}{C = {{YX}^{T}\left( {XX}^{T} \right)}^{- 1}} & (3) \\{Q = {\frac{1}{D}\left( {Y - {CX}} \right)\left( {Y - {CX}} \right)^{T}}} & (4)\end{matrix}$

where X and Y are the matrices formed by tiling the D total data pointsx_(t) and y_(t). For the Kalman filter, we also assume that the dynamicsof observed arm kinematics match the desired neural cursor kinematics,and so the parameters of the dynamics or trajectory model, A and W, arefit from observed arm kinematics:

$\begin{matrix}{A = {X_{2}{X_{1}^{T}\left( {X_{1}X_{1}^{T}} \right)}^{- 1}}} & (5) \\{W = {\frac{1}{D - 1}\left( {X_{2} - {AX}_{1}} \right)\left( {X_{2} - {AX}_{1}} \right)^{T}}} & (6)\end{matrix}$

X₁ is all columns of X except for the last column and X₂ is all columnsof X except for the first column, introducing a one time-step shiftbetween the two matrices.

In practice we constrain the form of the A and W matrices to obey simplephysical kinematics, integrated velocity perfectly explains position:

$\begin{matrix}{A = \begin{pmatrix}1 & 0 & {t} & 0 & 0 \\0 & 1 & 0 & {t} & 0 \\0 & 0 & a_{{vel}_{horiz},{vel}_{horiz}} & a_{{vel}_{horiz},{vel}_{vert}} & 0 \\0 & 0 & a_{{vel}_{vert},{vel}_{horiz}} & a_{{vel}_{vert},{vel}_{vert}} & 0 \\0 & 0 & 0 & 0 & 1\end{pmatrix}} & (7)\end{matrix}$

After fitting with either set of kinematics, a_(vel) _(vert) _(,vel)_(horiz) and a_(vel) _(horiz) _(,vel) _(vert) are typically close to 0and a_(vel) _(horiz) _(,vel) _(horiz) and a_(vel) _(vert) _(,vel)_(vert) are less than 1. The resulting model introduces damped velocitydynamics. Therefore, given no neural measurements, we expect a cursor inmotion to smoothly slow down. We also constrain the W matrix, so thatfor the dynamics model, integrated velocity fully explains position:

$\begin{matrix}{W = \begin{pmatrix}0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & w_{{vel}_{horiz},{vel}_{horiz}} & w_{{vel}_{horiz},{vel}_{vert}} & 0 \\0 & 0 & w_{{vel}_{vert},{vel}_{horiz}} & w_{{vel}_{vert},{vel}_{vert}} & 0 \\0 & 0 & 0 & 0 & 0\end{pmatrix}} & (8)\end{matrix}$

If we fit the full C matrix, the resulting filter is a position/velocityKalman filter (neural firing simultaneously describes position andvelocity). If we constrain the position terms to be 0, the resultingfilter is a velocity only Kalman filter (neural firing describes onlyvelocity).

Kalman Filter Design

Imagine a hypothetical prosthetic that decodes from a single neuron.This neuron fires more vigorously when leftward velocities areinstructed and also happens to fire more when the cursor is positionedon the left side of the workspace. If our decoder translates this firingrate into velocities, without any knowledge of positional effects, everytime the cursor enters the left side of the screen positive feedbackwill accelerate the cursor to the left. Positive feedback resultsbecause the firing rate increase due to leftward position is mapped toleftward velocity by the decoder, thereby driving the cursor faster tothe left than the user intends. By accounting for position, some of theincreased firing can be explained by the position of the cursor, andthis feedback effect can be mitigated.

However, the way in which the position/velocity Kalman filter models therelationship between position and velocity leads to undesired highfrequency jitter in the cursor position. The dynamics model described inthe Kalman filter section supra is physically based, acting like anobject moving in a gravity free 2-D space with damped velocity, so wemay expect cursor position to evolve smoothly. However, the Kalmanfilter translates velocity uncertainty into position uncertainty atsubsequent time-steps. To understand why this occurs, we examine how thefilter is run online. At time t we have a previous estimate of thekinematic state, {circumflex over (x)}_(t−1) and a new neuralobservation, y_(t). The first step in each filter iteration is to applythe dynamics model, estimating x_(t)=[p_(t), v_(t)] with all neuralobservations up to time t−1. This is the a priori estimate of x_(t):

{circumflex over (x)}_(t|t−1)=A{circumflex over (x)}_(t−1)  (9)

The model also estimates the a priori covariance (or uncertainty) of{circumflex over (x)}_(t|t−1):

P _(t|t−1) =AP _(t−1) A ^(T) +W  (10)

We know that the W adds no uncertainty to a priori position, given itsstructure as defined in Equation 8. However, AP_(t−1)A^(T) translatesprevious velocity uncertainty into current position uncertainty. Thismakes sense: if we do not know the previous velocity with certainty, wedo not know the integrated velocity with certainty and so our positionestimate may have error. Thus, in practice, there is uncertainty in thea priori estimate of every kinematic variable. To see how thisuncertainty in position translates to jitter in the decode, we cancontinue to step through the algorithm. Once we update the a prioriestimate, we must incorporate the information acquired from the neuralobservation. The model has an expected neural output given {circumflexover (x)}_(t|t−1), and this output may not match y_(t). This error isthe measurement residual, {tilde over (y)}_(t), and also has acorresponding covariance (or uncertainty) estimate, S_(t):

{tilde over (y)} _(t) =y _(t) −C{circumflex over (x)} _(t|t−1)  (11)

S _(t) =CP _(t|t−1) C ^(T) +Q  (12)

If this residual is nonzero (which is almost always true in practice),then the measurement and {circumflex over (x)}_(t|t−1) do not agree andwe must decide how much weight this observation residual has relative to{circumflex over (x)}_(t|t−1). This weight is based on how muchuncertainty is present in the kinematics suggested by the a prioriestimate of x_(t) versus the kinematics suggested by {tilde over(y)}_(t):

K_(t)=P_(t|t−1)C^(T)S_(t) ⁻¹  (13)

Finally, we can use K_(t), called the Kalman gain, to find the estimateof x_(t) that incorporates all of the neural observations up to time t,this is called the a posteriori estimate. We calculate the a posterioriestimate for x_(t) and its covariance:

{circumflex over (x)} _(t) ={circumflex over (x)} _(t|t−1) +K _(t){tilde over (y)} _(t)  (14)

P _(t)=(I−K _(t) C)P _(t|t−1)  (15)

The Kalman gain transforms the measurement residual into the kinematicspace. Since a priori estimates of both position and velocity kinematicshave uncertainty and neural measurements have information about positionand velocity, the Kalman gain will translate neural measurements intoupdated a posteriori estimates of both position and velocity. Foroffline trajectory reconstruction, this makes sense, as this coupledposition/velocity uncertainty exists throughout time. However, theseassumptions breakdown in the online setting, and substantially limitperformance.

We must distinguish online and offline use of the Kalman filter. In theonline setting, the user is presented with the a posteriori estimate ofcursor kinematics at every time-step. If we believe that the user seesand internalizes the presentation of the cursor on the screen at eachtime-step, then the way in which we model a posteriori covariance nolonger makes sense, as the user accepts the presented position as thecurrent position state. By presenting the decode to the user, we arecreating a causal intervention, that explicitly sets the value of thekinematic variable. This operation is defined by probability theory andis well described by the causal calculus.

As a first step to modify the filter to incorporate this feedback, wepresume that the user internalizes the filter's estimate of cursorposition, {circumflex over (p)}_(t−1), with complete certainty by timet. Accordingly, p_(t) is explicitly set to A_(p){circumflex over(p)}_(t−1), where A_(p) is the upper-left block of the state matrix A.In other words, since the user knows the previous cursor position{circumflex over (p)}_(t−1) via feedback, this forward model is exact atp_(t)=A_(p){circumflex over (p)}_(t−1). This is shown graphically inFIG. 8, where the intervened variable is shown by 810 (adding anothercircle is standard notation for causal interventions). Note also thatthe arrows coming into p_(t) have been removed, to indicate that p_(t)has been externally set and uncertainty is not propagated.

The result of this intervention is to remove uncertainty in p_(t). Allparameter fitting methods described in previous sections remainunchanged. To implement this position feedback filter, only a smallchange in the online operation of the standard filter is necessary. Allsteps outlined above are the same except for Equation 10. Previously, wehad:

$\begin{matrix}{P_{t|{t - 1}} = {{{{AP}_{t}A^{T}} + {W\mspace{14mu} {where}\mspace{14mu} P_{t|{t - 1}}}} = \begin{bmatrix}P_{p,p} & P_{p,v} \\P_{v,p} & P_{v,v}\end{bmatrix}}} & (16)\end{matrix}$

where each block of the matrix P_(t|t−1) represents the uncertaintypropagated from previous kinematic estimates (position to position,position to velocity, and so on). Since we have intervened and set p_(t)with feedback, this matrix becomes:

$\begin{matrix}{{\overset{\_}{P}}_{t|{t - 1}} = {\begin{bmatrix}0 & 0 \\0 & P_{v,v}\end{bmatrix}.}} & (17)\end{matrix}$

Otherwise, this filter is run in the same manner as the standard Kalmanfilter.

Brain-Machine Interface Utilizing Interventions to Emphasize Aspects ofNeural Variance and Decode Speed and Angle

BMI systems utilize neural recordings to drive prosthetic devices. Theperformance of BMI systems is important to their clinical viability.Specifically, the highest performing BMIs measure action potentials(spikes) from neurons in motor cortical regions. These systems rely onneural observations from electrode arrays implanted into motor corticalregions. However, the quality of the recordings on the electrodes isknown to degrade over time, which in turn degrades the quality of theBMI. In this embodiment of the invention we describe a technique toimprove BMI performance after it has degraded significantly due to lossof neural recordings (FIG. 9).

Experimental Approach

In the spirit of addressing the objective of this invention, we modeled(1) a poorly performing decoder, and then subsequently (2) a decodertrained after 25% of available channels (neural recordings) are randomlylost. In one example, we showed that after having lost 25% of channelsso that state-of-the-art decoders perform very poorly, the invention(called the intervention PCA algorithm) can recover, and even exceed theperformance of state-of-the-art decoders when 100% of the channels wereavailable. Specifically, we compared the intervention PCA algorithm tothe highest-reported performing algorithm for 2D cursor control, theReFIT-KF (Gilja et al., Nature Neurosci, 2012).

Evaluation Task

Performance was evaluated with the grid task, where the subject canselect one of several targets in a grid. This task allows us to quantifyperformance in terms of the amount of information conveyed in a giventime, or the achieved communication rate (in bits per second). Moreover,the task is clinically relevant, since the grid mimics a keyboard. Weused the optimal keyboard layout found via optimization, which is a gridof 25 targets with a selection time of 450 ms.

Quality of Observations

To bias ourselves in a regime where decoder performance is poor, werandomly selected N=48 channels from the array as reporting observationsof neural spiking This resulted in a poorly performing decoder. Then, weevaluated the decoder performance after randomly discarding 25% of thepreviously selected 48 channels, so that there remained N=36 channels.We evaluated the ReFIT-KF (state-of-the-art) when 100% of channels wereavailable and when 75% of channels were available. We evaluated theintervention PCA algorithm (invention) only when 75% of the channelswere available.

Choice of Algorithm

The intervention disclosed in this invention can be used with anydecoder in which we can easily control the parameters, which map neuralobservations into prosthetic kinematics. This includes the class oflinear BMIs (e.g. simple least-squares regression, Wiener filters,Kalman filters, etc.).

To evaluate what the baseline performance when 100% and 75% of channelswere available, we used the ReFIT-KF algorithm (Gilja et al, 2012). TheReFIT-KF algorithm is the highest-performing algorithm for 2D cursorcontrol available in the literature, by almost a factor of 2×. The newalgorithm of this invention to recover performance when 75% of thechannels are available can be called “intervention PCA-ReFIT KF,” or“intervention PCA” for brevity. For the intervention PCA decoder, weused as observations the low-dimensional projections (found via PCA) ofthe neural data. This means that the activity of N channels is projectedinto a D dimensional space, where D<N, and treated as the neuralobservations of the BMI. The details of the intervention PCA algorithmare discussed in section 2. Inspired from the intervention PCA approach,we also discuss how a decoder can be modified to output, in addition tothe position and velocity of the prosthesis, the speed and direction ofthe prosthesis.

Intervention PCA

The intervention PCA algorithm leverages scientific knowledge of thelow-dimensional nature of neural data. For the intervention PCAalgorithm, we modified the ReFIT-KF algorithm by making the observationsof the Kalman filter the principal components of the data (as opposed tobinned spike counts) and performing an intervention on the noise models,which is discussed in more detail below.

Projection to Reproduce Original Trajectories

When 100% of the channels were available, we used dimensionalityreduction (via PCA) to compute the low-dimensional “canonical” neuraltrajectories. These trajectories represent the projection of all theneural activity into a lower-dimensional space. When only 75% of theavailable channels remained, we performed a regression so that theremaining 75% of neural data attempted to reproduce the canonicaltrajectories (remembered from when 100% of the neural data wasavailable). Therefore, we learned a matrix which mapped the neural datainto a low-dimensional space (albeit noisily) so that the neuraltrajectories (which are the observations of the BMI) approximate theneural trajectories when 100% of the channels were available. This wasperformed by least-squares regression.

Move-stop Variance is Largely Contained in the Top Principal Componentsof the Data

Whether the subject is moving or stopping the prosthetic device is alarge variance event, evidenced by a general rise in firing rates acrossmany channels. Because this event (capturing the speed of theprosthesis) is of large variance, and present in a significantproportion of channels, much of this variance is present in the topprincipal components (PCs) of the neural data since they reflectdimensions capturing the greatest amount of neural variance.

In linear BMI systems, this means that PC₁ typically contributes a largevelocity to the decoder. However, due to the limitations of lineartechniques used for BMIs, if one is decoding velocity, then PC₁contributes a large magnitude speed in only one direction. In the caseof decoders with lower performance, this can severely bias the decoder,if there is not much directional tuning in PC₁. This is because a largeprojection of the data in PC₁ will drive the decoder in a certaindirection with significant magnitude, even though it has littledirectional tuning

The Kalman Filter Weighs Observations in Accordance to the ObservationProcess Noise

Although this technique could be applied to other types of linearfilters, we will discuss its particular application to Kalman filters(and will briefly mention an example of how to achieve a similar effectin other linear techniques). We specifically address Kalman filtersbecause they are the highest-performance algorithms reported inliterature for 2D cursor control.

State-of-the-art BMIs are based on velocity Kalman filters (VKFs), whichdecode an underlying state x_(k) (a vector denoting the velocity of thecursor), and the observations y_(k) (a vector denoting the neuralobservations). These Kalman filters operate on a linear dynamical system(LDS) of the form:

x _(k+1) =Ax _(k) +w _(k)

y _(k) =Cx _(k) +q _(k)

where w_(k)˜N(0,W) and q_(k)˜N(0,Q). For intervention PCA, where the(y_(k)) are the principal components of the neural data, (y_(k))_(i)refers to PC_(i).

The observation noise covariance, Q, models the covariance in theobservation process, y_(k)−Cx_(k). Hence, if Q₁₁ is relatively large, itheuristically indicates that estimating the first dimension of theobservation, (y_(k))₁ with the underlying state, (C_(1,:))x_(k), isnoisy (where C_(1,:) is the first row of C). Therefore, the firstdimension of the observation is relatively unreliable. The interventionwill decrease the reliability of the observation of PC₁ which decreasesits contribution to the decoder. With the VKF, this is possible byincreasing Q₁₁. In other linear techniques, e.g. if there is a singlematrix mapping L which calculates x_(k)=Ly_(k), it is possible todownweight the observation (y_(k))₁, by scaling the coefficients in thefirst column of L by a constant between 0 and 1.

The Intervention PCA Decoder can Increase the Performance of PoorlyPerforming Decoders

The goal of the intervention PCA decoder is to mitigate degradingdecoder performance, specifically by ameliorating biases that typicallycripple decoders. In our particular experiments, intervention PCAaddresses this by increasing the observation noise of PC₁ in a Kalmanfilter (although we allow the possibility of doing this for otherprincipal components, e.g. PC₂). This effectively downweights thecontribution of PC₁ to the decoder. As previously mentioned, the mainmotivation for downweighting PC₁ is because PC₁, through the KF,contributes a large speed in a certain direction by virtue of being alarge variance event associated with moving and stopping of the cursor.Hence, even if PC₁ has little directional tuning, due to the Kalmanfilter modeling, it contributes a massive speed in a single direction bythe modeling of the Kalman filter. Thus, by downweighting thecontribution of PC₁, we ameliorate any major bias that may arise fromsuch modeling assumptions. We reiterate that this technique is notlimited to Kalman filters.

It should be noted that the intervention PCA decoder is not merelydownweighting the observations relative to the kinematics model of aVKF, x_(k+1)=Ax_(k)+w_(k). That is, it is not merely increasing themomentum of the cursor. In fact, the amount of weight given toobservations relative to the kinematic model is approximately equalbetween intervention PCA and ReFIT-KF. Instead, what the interventionPCA decoder is doing is differentially weighing aspects of neuralvariance. If PC₁ does not have large directional variance, but largespeed (move/stop) variance, then downweighting PC₁ would allow otherdimensions, which may have more directional tuning, to be more prominentin the decoder. Hence, by using PCA, we are able to do a globaldownweighting (affecting all channels), varying for each channel, asthey contribute to the dimension of interest. In this sense, thelow-dimensional projection can be viewed as approximately partitioningaspects of neural variance into different dimensions; we then perform anintervention to downweight certain dimensions. It should be noted thatother dimensionality reduction techniques which separate aspects ofneural variance could be used, so that the technique is not limited toPCA. Although PCA is not guaranteed to partition the speed (move-stop)vs directional variance amongst different dimensions, we have causallyshown that PC₁ contains a sufficient amount of information to determinewhether the subject is moving or stopping, and is the largestcontributor to velocity. Hence, there is approximate varianceorganization resulting through PCA, which allows us to differentiallyattenuate a dimension with less directional tuning (and more to do withspeed variance) relative to other dimensions. This increases theperformance of poor decoders.

Intervention PCA Results

First, we showed (FIG. 10) that the performance of ReFIT-KF degradesafter dropping 25% of the channels randomly. We evaluated theperformance of ReFIT-KF when 100% of channels were available. In thefigures below, where red denotes the performance of the state-of-the-artReFIT-KF, each small dot denotes the bitrate of an experimental blockwhere the subject used a decoder to select targets on a grid, while thelarge dot denotes the mean bitrate across all experimental blocks. When100% of the channels were available, the variance in performance was dueto variance in subject performance within a day, and then across days.Then, we evaluated the performance of ReFIT-KF when x=75% of channelswere available, to show that the performance degrades with more channelloss. For each experimental block (small dot), we randomly selected 25%of the channels and subsequently dropped them, so that each block has adifferent set of channels dropped. We note that dropping 25% of thechannels leads to approximately 60% of the performance on the grid task(compared to when 100% of the channels were available).

Next we show (FIG. 11) the performance of intervention PCA (1110) forthe same experimental blocks where 25% of the channels were randomlylost. We note that intervention PCA recovers performance, and indeedsurpasses the performance of ReFIT-KF at 100% of the channels.

Decoding Speed and Angle

Inspired by the intervention PCA results, we also evaluated thefeasibility of building a decoder that regresses the neural data againstthe desired speed and angle of the prosthesis (as opposed to theCartesian velocity and position) and outputs these variables to akinematic space. Thus, the dimensions of the output kinematic space(i.e., speed and angle) are independently and separately estimated suchthat they are not derived from Cartesian velocities. Note that prioralgorithms would necessarily have to combine the independent dimensionsestimated by the decoder to arrive at speed and angle, since they do notestimate these variables separately and independently.

Thus, we separately and independently decoded the speed and angle of theprosthesis rather than the Cartesian velocity of the prosthesis. In oneembodiment, we separately and independently decoded the speed by makingit the state of the Kalman filter (rather than the Cartesian velocity).After decoding a speed variable, s, we allowed the speed to betransformed according to a non-linear function, s^(n)=f(s). In oneembodiment, this non-linear function was a piece-wise linear functionthat amplified higher speeds. To decode the angle, θ, of the prosthesis,in two embodiments, we decoded both the sin(θ) and cos(θ) of theprosthetic device using either a Wiener filter or a Kalman filter. Thenthe direction of movement could be calculated as θ=arctan(sin(θ)/cos(θ))or the x- and y-velocities could be calculated as s^(n) cos(θ) and s^(n)sin(θ) respectively.

REFERENCE

-   (1) V. Gilja, P. Nuyujukian, C. A. Chestek, J. P. Cunningham, B. M.    Yu, J. M. Fan, M. M. Churchland, M. T. Kaufman, J. C. Kao, S. I.    Ryu, and K. V. Shenoy, “A high-performance neural prosthesis enabled    by control algorithm design,” Nature Neuroscience, vol. 15, pp.    7-10, Nov. 2012.

What is claimed is:
 1. A brain machine interface, comprising: (a) analgorithm executable by a computer for separating neural signals into:(i) neural signals related to a speed of a prosthetic device and (ii)neural signals related to an angle of movement of the prosthetic device;and b) a decoder algorithm executable by the computer for mapping thespeed-related neural signals and mapping the angle-related neuralsignals to a kinematic space or the prosthetic device.
 2. A brainmachine interface, comprising a decoder algorithm executable by acomputer for regressing neural signals against an angle and magnitude ofa velocity and outputting to a kinematic space or a prosthetic devicesuch that the dimensions of the output kinematic space or the prostheticdevice are independently and separately estimated and not derived by acombination of the output of the decoder algorithm.
 3. A brain machineinterface, comprising a decoder algorithm executable by a computer forincorporating previously known information regarding neural signals thatare no longer present due to the loss of neural recording, neuron celldeath, or noise.